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We present a continuous-time Monte Carlo impurity solver for multiorbital impurity models which 
combines a strong-coupling hybridization expansion and a weak-coupling expansion in the Hund’s 
coupling parameter J. This double-expansion approach allows to treat the dominant density-density 
interactions U within the efficient segment representation. We test the approach for a two-orbital 
model with static interactions, and then explain how the double-expansion allows to simulate models 
with frequency dependent U{ijj) and J(cu). The method is used to investigate spin state transitions 
in a toy model for fullerides, with repulsive bare J but attractive screened J. 

PACS numbers: 71.10.Fd 


I. INTRODUCTION 

Strongly correlated materials exhibit a range of in¬ 
teresting properties, such as unusually large suscepti¬ 
bilities or high-temperature superconductivity. A the¬ 
oretical investigation of this class of materials is possi¬ 
ble within the framework of dynamical mean-field the¬ 
ory (DMFT)ji either at the simple model level or in 
combination with input from density-functional based ab 
initio calculations!^ Because of the interest in cuprate 
high-temperature superconductors, and also because of 
algorithmic limitations, much effort has in the past 
been devoted to the study of the single-band Hubbard 
model. The discovery of aromatic superconductors,— 
iron pnictides^ and recent studies emphasizing the multi¬ 
band character of cuprates^ have however shifted some 
of the attention to correlated multi-orbital systems. In 
these multi-orbital systems, the Hund’s coupling param¬ 
eter plays an essential role and leads to new types of 
correlation phenomena, such as bad-metal behavior due 
to local-moment formatioii)^^— magnetism and orbital 
ordering^^— spin-state transitionsorbital-selective 
Mott transitioned^ nontrivial spatial correlations 
and unconventional superconductivity—^ 

Dynamical mean field simulations of generic multi¬ 
band models have become possible thanks to the develop¬ 
ment of strong-coupling (hybridization expansion) impu¬ 
rity solvers— S The matrijs^i or Krylov-implementatione 
of this impurity solver can handle arbitrary interactions 
among the orbitals, but the computational effort scales 
exponentially with the number of orbitals. A far more 
efficient simulation is possible within the so-called seg¬ 
ment formalism^dS if the interactions are restricted to the 
density-density component of the full Coulomb matrix, 
and this approximation is still often made in simulations 
of transition metal and actinide compounds. In most 
cases these density-density terms give the dominant con¬ 
tribution to the interaction energy, so that it may be 
advantageous to consider the spin-flip, pair-hopping, and 
correlated hopping terms as a perturbation in an expan¬ 


sion that treats the density-density components exactly. 
In this paper, we explore such a double-expansion im¬ 
purity solver, which stochastically samples a diagram¬ 
matic expansion of the impurity partition function in 
powers of the hybridization function and the interaction 
terms which are not of density-density type. Such an 
impurity solver trades the exponential scaling of the ma¬ 
trix/Krylov approach with an additional weak-coupling 
type expansion, and (for more than two orbitals) a po¬ 
tential sign problem. It should be efficient in the case of 
a small number of orbitals and not too large Hund’s cou¬ 
pling. Here, we implement and test the double-expansion 
solver for a two-orbital model with rotationally invariant 
interaction. 

Besides potential efficiency gains, a second important 
reason for exploring the double-expansion approach is 
that such a solver enables the simulation of certain types 
of problems which cannot be solved using the established 
hybridization expansion methods. A relevant example 
is models with a dynamically screened J. The low- 
energy effective models solved in DMFT simulations of 
correlated materials can be obtained from a down-folding 
procedure in which the bands outside some energy win¬ 
dow around the Fermi level are integrated outi^ This 
procedure leads to a dynamically screened interaction. 
For example, in transition metals and their compounds, 
the screening of charge fluctuations results in density- 
density interactions which range from a bare value of 
typically about 20 eV to a screened value of only a few 
eV. Highly efficient algorithms exist to treat this type of 
screening— The screening of the Hund’s coupling pa¬ 
rameter is usually much weaker, so that the bare and 
screened J typically differ by less than 20%i^ In all 
the DMFT simulations to date, the Hund’s coupling has 
thus been treated as frequency-independent. However, 
given the sensitivity of multi-orbital phase diagrams on 
the Hund’s coupling parameter— it is desirable to 
develop a method which extends the efficient technique 
of Ref. [2^ to models with a frequency-dependent J{uj). 
There are also materials in which the dynamical screening 
of J plays a crucial role. In alkali-doped fullerenes^i^S 
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the observed superconductivity is believed to arise from 
an overscreened as a result of Jahn-Teller screen¬ 

ing J(uj) turns negative at some low frequency and hence 
favors low-spin states. The double-expansion solver al¬ 
lows simulations of such dynamically screened multi¬ 
orbital systems. 

The structure of the paper is as follows. In Section m 
we describe the double expansion method for the case of 
a two-orbital model with rotationally invariant interac¬ 
tions and explain how this method can be used to treat 
models with dynamically screened U and J. Section IIIII 
shows some test results and information on the average 
perturbation orders. Simulation results for a model with 
static U and dynamical J(w) are presented in Section HVl 
and a brief summary and outlook is given in Section |Vl 


number k. The bath energies and the hybridization pa¬ 
rameters Vk^a,a define the hybridization function 
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In the case of a semi-circular density of states with band¬ 
width At a, the DMFT self-consistency condition provides 
a simple relation between the hybridization functions and 
impurity Green’s functions Ga,cr'^ 
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We will consider an orbital-independent semi-circular 
density of states with ta = t {a = 1,2) and use t as 
the unit of energy. 


II. MODEL AND METHOD 


A. Two-orbital model 


B. Double expansion for static interactions 


As a simple, but nontrivial example, we consider 
a two-orbital model with rotationally invariant Slater- 
Kanamori interactions and a possible crystal field split¬ 
ting. DMFT replaces the lattice problem by the self- 
consistent solution of a two-orbital quantum impurity 
model with Hamiltonian 

H = 'Hdens + Hsf + nlf -I- "Uph + Tfpjj 


+7fbath + Tfhyb + 

(1) 

The density-density, spin-flip, pair-hopping, bath and hy¬ 
bridization parts of the Hamiltonian are given by 

^dens — ^ ^ ^ ^ 

a,O' <T 

- n2,cr) 

Q! CT 

^1,(T^2, —(T 

“f" ^ y^ 

(2) 

^Sf *^^1,4, ^2,4,^1,'I', 
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(4) 

^bath — ^ ^ 

k,06,a 

(5) 

k,06,a 

(6) 

where a = 1,2 is the orbital index, a =td (or ±1) 
the spin index, A the crystal field splitting, U the intra- 


orbital interaction, and J the coefficient of the Hund cou¬ 
pling. We choose the interorbital interaction U' = U — 2J 
for rotational invariance and denote the impurity cre¬ 
ation operators by cj, g. and the density operators by 
na,(j = cl^g.Ca,a- The bath levels, with creation opera¬ 
tors Q, o- energy Ck are parametrized by a quantum 


We solve the impurity model using the continuous¬ 
time Monte Carlo techniquei^ In the double-expansion 
approach, this continuous-time method is based on a si¬ 
multaneous expansion of the partition function in the hy¬ 
bridization terms and the interaction terms which are not 
of density-density type (in the model considered here, the 
spin flip and pair hopping terms). To derive the formal¬ 
ism, we switch to an interaction representation in which 
the time-evolution of operators is given by 'Hdens + bath 
and write the partition function of the impurity model 
as 


Z = Tt. 


Trf, e ^(’^dena+Wbath) y 

+ + 'Hphir) + + 'Hhyb('r) -f j 


(9) 


The next step is to expand the time-ordered exponential 
in powers of the spin-flip, pair-hopping and hybridization 
terms. Since there is then no coupling between the impu¬ 
rity and the bath anymore in the time-evolution (given 
by 'Hdens + Tfbath), the trace over the bath states can be 
computed analyticallyi^i^i This leads to the expression 



where Z^ath = Tr^e and the weight of a con¬ 

figuration consisting of 2n/i hybridization events (rih = 
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FIG. 1: Illustration of an order Uh = 12, rist = 2 configuration 
for the two orbital model. Empty (full) circles represent an¬ 
nihilation (creation) operators. Vertical dashed lines indicate 
spin-flip events. 
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FIG. 2: Illustration of an order nu = 7, Uph = 2 configuration 
for the two orbital model. Empty (full) circles represent an¬ 
nihilation (creation) operators. Vertical dashed lines indicate 
pair-hopping events. 


2nsf Spin-flip events and 2nph pair hopping 
events is given by 


w = Tlr 


3 /3’Hdens 
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C. Monte Carlo sampling 


The expression for the trace in Eq. (HI) shows that 
a given configuration C consists of 2nh hybridization 
events, 2nsf spin-flip events, and 2nph pair-hopping 
events. We can generate all possible configurations using 
local updates which insert or remove pairs of hybridiza¬ 
tion, spin-flip or pair-hopping events. These updates 
must satisfy the detailed-balance condition w{C)p{C —f 
C) = w{C')p{C' —>■ C), where p{C C) is the tran¬ 
sition probability from configuration C to configuration 
C. We split this transition probability into a proposal 
probability and an acceptance probability, p{C —>■ C) = 
pP™P(C C")p^“(C' C) and define the ratio of ac¬ 
ceptance probabilities 


where the time evolution of operators is now given by 
^dens; S — ^ and 

is a na,cr x na^a matrix of hybridization functions, 
with elements M“^(i, j) = AA^hi “ 

The trace vanishes unless there is an equal number 
of impurity creation and annihilation operators for each 
flavor. This requirement implies that for each T^hyb, we 
must have a corresponding and similarly for the 

spin-flip and pair-hopping terms. The expansion in the 
spin-flip and pair-hopping terms hence does not lead to 
a sign problem in this two-orbital case with static in¬ 
teractions. Furthermore, because the time-evolution op¬ 
erator is diagonal in the occupation number basis, we 
can use the segment representations^ to graphically rep¬ 
resent all the non-vanishing contributions to the trace 
(see Figs. [Hand[2|). In this representation, each segment 
marks a time-interval in which the impurity is occupied 
by an electron of a given flavor (a, a). We denote such a 
segment configuration by C and sample the space of all 
configurations using the Metropolis algorithm. 


B(C ^ C') = ^ ^ C) w{C') 

1 > paccj'j^/ ^ C') pP™P(C' ^ C) W{C) ■ 

( 12 ) 

In the Metropolis scheme, the update is accepted with 
probability min[l, R], 

For the hybridization events, the sampling procedure is 
exactly the same as detailed in Ref. [l^- For the insertion, 
we try to place a creation operator at a randomly chosen 
time A Ih® imaginary-time interval. If it falls on a seg¬ 
ment, the move is rejected. Otherwise, we compute the 
length /max of the interval to the next creation operator 
(which may be associated with a spin-flip or pair-hopping 
event) and choose the time for the annihilation oper¬ 
ator randomly in this interval. (By next operator we 
always mean the neighboring operator in the direction of 
increasing imaginary time, taking into account periodic 
boundary conditions.) The distance between the inserted 
operators will be denoted by 1. For the removal, we ran¬ 
domly pick one of the creation operators, and propose to 
remove the segment attached to this operator, provided 
it is not cut by a spin-flip or pair-hopping term. The 
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corresponding acceptance ratio reads 


Rhyh{n^a,(T ^ “t“ 1) — 
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(13) 


where = ±A for a = 1,2 and I overlap denotes the 
total length of the overlap between the inserted segment 
and the segments associated with the /3, a'-state. 

In addition to the hybridization expansion algorithm 
updates, we sample the spin-flip and pair-hopping terms 
in C using the following updates: 

1. insertion and removal of 5'(rs)5''l'(T(), Tg > r(. 


2. insertion and removal of (t( ) S'(ts ), > Ts, 

3. insertion and removal of P{Tp)P^{Tp), Tp > r^, 

4. insertion and removal of P^Tp)P{Tp), > Tp. 

The insertion of spin-flip and pair-hopping events is only 
possible if the impurity is in the appropriate local state, 
namely the |i, t)j Itil); |0:ti) and It)-, 0) in the order 
given above. 

Suppose we want to insert a spin-flip operator 
S{Ts)S^{Tg) with Ts > Tg. First, we generate a random 
imaginary time r' and check whether the local state at 
Tg is the appropriate one, namely t) in this particular 
case. If it is not the correct local state, the move is re¬ 
jected. If the insertion is possible, we compute the length 
^max from Tg to the next operator and choose Tg randomly 
within this interval. In the inverse procedure, we remove 
an SS^ operator. To do so, we randomly select an 
operator among the Ugf operators and remove it to¬ 
gether with the S operator next to it (in the direction of 
increasing time), provided that there are no segments or 
pair-hopping operators in between. The same strategy is 
used for the insertion/removal of the 5'^' (t( ) S'(ts ) opera¬ 
tor, and for the sampling of the pair-hopping operators. 

If we denote the length of the interval between and 
Ts by Z, the acceptance ratios for the spin-flip and pair¬ 
hopping events become 

Rssi{.nsi^ns{+1) = (14) 

risf + 1 

Ugi + I 

i?PPt (riph -)■ riph + I) = (16) 

Tlph “T i 

flptp(nph-» riph + 1) = (17) 

(Here, we assume that the operator on the right has the 
smaller time argument.) The acceptance ratio for the 
two spin-flip events is the same because the intra-orbital 
interaction is spin symmetric, and the occupation of the 
orbitals does not change. The asymmetry in the pair¬ 
hopping case comes from the crystal field splitting, which 
favors the occupation of one of the orbitals. 


D. Retarded interactions 


As mentioned in the introduction, realistic low-energy 
models of correlated materials involve retarded interac¬ 
tions. In the case of the two-orbital model considered 
here, the DMFT impurity action thus takes the follow¬ 
ing general form 


5'efl — So + Sint + <5'] 


hyb, 


(18) 


where So contains the chemical potential and crystal field 
terms, i7hyb the hybridization functions and Sint is given 

by 


'S'int = ;r^ / dr / dr'na.<T(r)Wf;^ (r - T')n;3,o./(r') 


a,(7 
/3,ct' 


oX! / dT / dT'J(r-T')fc^(T)A2)(T') 


/o ^0 


+ A2i(r)X]?(r') +Xi2(r)X]?(r') (r)A^)(r') 

(19) 


where = c|_^C 2 ,ct and = cj 

The time-dependent interactions have an instanta¬ 
neous component corresponding to the bare interaction 
and an (attractive) retarded component describing the ef¬ 
fect of screening: Jir') = Jbare<5(r) -I- Jret(r) and U{t') = 
t7bare<5(T)-I-t/ret(r). Note that the Jbare<5(T) Contribution 
in the second term yields not only the instantaneous spin- 
flip and pair-hopping terms 'Hsf + 'H\e + 'Hph + (with 
J = >7bare)j but in addition also a same-spin inter-orbital 
density-density interaction — Aaare X)a<,s a-''^a.o-^/3,0' ^-'^d 
a chemical potential term ^ Aaare o-^a,cr- Therefore, 
in the rotationally invariant case, the density-density in¬ 
teractions in the first term are chosen as 

K%.{t)=U{t\ (20) 

= UlllLr) = U{t) - 2J{t). (21) 

The retarded interactions C/ret('r), Jret(T) arise from 
some electron-boson coupling term of the form 

= E E + bl), (22) 

2^ a,Q;',(7 

where fct denotes the creation (annihilation) oper¬ 
ator for the izth bosonic degree of freedom. The bo¬ 
son one-body part is given by H^re, the 

bosonic degrees of freedom represent plasmons, bosonic 
modes corresponding to single-particle excitations, and 
phonons. The former two are responsible for the dynam¬ 
ical screening processes which are taken into account in 
the down-folding procedure. The phonons further reduce 
the resulting interactions values. 
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In terms of and Wj/, f/ret(w) are Jret(w) can be 

written as 


TT ( \ 

C/ret(w)-2^ 


2^ ,.a - 


(23) 


t(w) - ^ , ,2 _ ,,,2 


2(A2_]^)^ _ ^Aj 2A2 _i 

m 2 — , .2 /,,2 — /,,2 


(24) 


Nowadays, ab initio estimates of these retarded interac¬ 
tions can be obtained by the constrained random phase 
approximation (cRPA)^ for the electronic screening con¬ 
tribution, and by the constrained density-functional per¬ 
turbation theory (cDFPT) ^^i^^ for the phonon contribu¬ 
tion. In reality, there is a continuum of electronic screen¬ 
ing frequencies, so that the sum in Eqs. (1^ and (IMl) 
becomes an integral over frequencies 

In most strongly correlated materials, the screened U 
and J remain positive, and the frequency dependence 
of J is rather weaki^ An interesting exception are the 
alkali-doped fullerides, where the Wannier functions are 
delocalized on a Ceo molecule. Here, the bare exchange 
interaction Jbare becomes small (Jbare ^ 0.1 eV, band¬ 
width ~ 0.5 eV)i^ The electronic screening contribu¬ 
tions reduce the static value of J to ~ 0.035 eV, which is 
still positive. The important low-energy screening con¬ 
tributions come from phonons, which yield an attrac¬ 
tion of ^ —0.05 eV, inverting the sign of the static ex¬ 
change interaction.— The resulting multiorbital Hamil¬ 
tonian with an overscreened J is of great interest, since 
a negative Jscr can induce nontrivial synergies between 
the correlated electrons and phonons, leading to an ex¬ 
otic s-wave superconductivity— r— However, the effect of 
the dynamical screening of J is still an open issue. The 
algorithm presented in this paper provides a basis for 
attacking this challenging problem. 

In solving the impurity model with the action SeS 
in Eq. CHI), we treat the density-density and spin- 
flip/pair-hopping terms in Eq. (1191) separately. It is 
therefore convenient to shift the Jha.TeS{T)S^^a-> contribu¬ 
tion to the density-density term and to write the ac¬ 
tion Sint with the interactions J and U replaced by 

J(t) = Jhnre5{T)5a,-cr' -f Jret(T) and W(r) = [/bare'5(r)- 

«/bare^(T)<lCT,cr'(l —^c(,/ 3 ) + Hret(T). In this formulation, the 
density-density interactions become 

K^a' (^) = C^bare^(T) + C/ret(T), (25) 

= ([/bare " 2Jbare)'5(T) 

+ ([/ret(r) - 2Jret(T)), (26) 

= ([/bare " 3Jbare)<5(T) 

+ ([/ret(r) - 2Jret(T)). (27) 


The half-filling condition, which is ^ 1/2 = §[/bare — 
fTbare in the model without dynamical screening, be¬ 
comes [l\j 2 — 2 ^scr 2Tscr with User ~ U{ui = 0) 

and Jscr = </(w = 0). In the limit Jret(7‘) = (Jscr — 


2T 

2i 

1T 

1i 



FIG. 3: (Color online) Illustration of an order rihyb = 4, Usi = 
1 configuration for a model with retarded density-density in¬ 
teractions. Dashed lines indicate interactions (v) con¬ 

necting all pairs of hybdrization, spin-flip and pair-hopping 
events (only the interactions involving the creation operator 
of the non-cut segment are shown). The sign SiSj associated 
with each dashed line is ±1 depending on the types of oper¬ 
ators (creation/annihilation) which it connects. 


•Aiare)<^(T) (high frequency screening) we recover the rota- 
tionally invariant 5'int with static interactions equal to the 
screened values. (Again, one has to take into account the 
additional density-density and chemical potential terms 
resulting from the and X^^X^^ operators.) 


1. Retarded density-density interactions 

The density-density contribution to S'int can be calcu¬ 
lated efficiently using the technique discussed in Refs. [13, 
The retarded part leads to an additional weight 
factor u>screen({Ti}) of the form 

Wscreen({ri}) = exp I ^ (^ - T, ) 

\ 2n>z>_?>l 

(28) 

where n = nt + Ugi Uph, the Ti are the times corre¬ 
sponding to segment start- or end-points (which can be 
the locations of hybridization, spin-flip or pair-hopping 
operators), and s = ±1 is a sign (-I-I for segment start- 
points and —1 for segment end-points). The function 
(r) is obtained from the twice integrated retarded 
interaction [/ret('r) and Jret('r)^ and thus is also orbital- 
and spin-dependent: 

K::'(.r) = Ku{r), (29) 

ir) = Kll, (r) = Ku{t) - 2Kj{t). (30) 

with = [/ret('r) and K'^{t) = Jret(T). Both Ku{t) 

and Kj{t) are defined in the range 0 < t < (3, are 
/3-periodic and symmetric around r = 5/2, and satisfy 
Ku,ji0+) = Ku,jW~) = 0. 

The structure of the configurations thus remains the 
same as in the simulations without retarded density- 
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density interactions (they consist of a collection of hy¬ 
bridization, spin-flip and pair-hopping events), the only 
difference is that now each pair of creation/annihilation 
operators (both hybridization events and segment start¬ 
er end-points corresponding to spin flip and pair hopping 
events) are linked by lines representing the “interaction” 
^a,a^p,a Fig.jS] In a local update, only 

the lines connected to the inserted or removed operators 
have to be considered. 

More explicitly, the acceptance ratio for a spin-flip (or 
pair hopping) insertion becomes 


^max 



Xt 




1 "^bare 


:xf 


0 

X x’ 

P 


-Rsst OC 

i?PPt OC ** 



? 

1 


(31) 

(32) 


where the sum over i runs over all operators associated 
with the new spin-flip (or pair-hopping) event, and the 
sum over j runs over all the other operators in the con¬ 
figuration pairs in the same orbital can be ignored). 


2. Retarded spin-flip and pair-hopping terms 

The sampling of the retarded spin-flip or pair-hopping 
operators is analogous to the algorithm discussed in 
Ref. Idol i.e. we expand the partition function in pow¬ 
ers of these terms and sample their contribution stochas¬ 
tically. In order to insert a retarded spin-flip operator, 
we replace an instantaneous spin-flip event by a retarded 
operator, with t ^ t' and cr yf a'. We 
describe the procedure here for retarded spin flips (the 
retarded pair hoppings are sampled in the same man¬ 
ner). First, we randomly select the type of operator {S 
or S'!) which we want to split. In the following, we as¬ 
sume it is an S operator. Next, we choose one of the nff 
instantaneous spin-flip events corresponding to S opera¬ 
tors. Suppose we select the S operator located at time r. 
This operator can be written as S{t) = X^'^{t)X^^{t). 
We randomly select either the X^ or X^^ operators for 
the proposed shift on the time axis (suppose it is X^'^). 
To fix the new location, we compute the distance Imax to 
the next operator in the forward-direction (taking peri¬ 
odic boundary conditions into account) and choose the 
time t' randomly within the interval of length /max- The 
proposed move is from the instantaneous spin-flip event 
with weight — JbareW|^(r)X.^^(r)c/r to the retarded spin 
flip with weight — Jret(T^ — t)X^^{T) dT^. In this 
new configuration, the hopping operators are connected 
by the interaction line Jret{T'—T) = J{t'—t) (see Fig.|l|). 

In the inverse move, we remove the retarded spin-flip 
event by randomly selecting one of the n® f retarded spin- 
flip pairs corresponding to S. We then randomly choose 
one of the operators and try to shift the other operator 
to its position on the time-axis. If there are other oper¬ 
ators between the retarded spin-flip operators, the move 
is rejected. With these procedures, the acceptance ratio 


FIG. 4: Retarded spin-flip insertion. An instantaneous spin- 
flip operator S{t) — X^^{t)X^^{t) is split into two separate 
hopping events X^^{t) and XI^(t') by randomly choosing the 
position of one of the operators (here Xj^) in the interval of 
length /max. 


for the retarded spin-flip insertion becomes 

Rxxirisf, nfsf nff — I, + 1) = 

_ w Jretid' "^) 

2-^ OL.a ^overlap_ 


‘rsf 


-hi 




bare 


(33) 


Because the number of instantaneous spin-flip events nff 
associated with S operators can now be different from 
the number associated with S'! operators, we have 
to keep track of these perturbation orders separately. In 
the instantaneous spin-flip updates, one then uses in 
Eq. (fT4l) and in Eq. (|T5|l . We also note that in the 
usual situation where Jbare > 0 and >7 (t) < 0 (0 < r < 
/3), Eq. (1331) implies that the Monte Carlo sampling for 
the simulation with retarded spin-flip and pair-hopping 
terms will suffer from a sign problem. 

If the screening frequency is high, so that Jret(T) 
approaches a (5-function, it is more efficient to absorb 
this factor into the proposal probability. More specifi¬ 
cally, we propose the time t' in the interval of length 
/max according to the probability distribution — Jret(T^ ~ 
t)/c/t'I Ji.et(T' — r)|. In this case, the ratio of 
acceptance probabilities becomes 


o c 


3. Kink updates 

Because of the retarded X„Xa terms in Eq. m , an er- 
godic sampling requires additional updates. One of these 
additional updates swaps the Jret-links of two randomly 
chosen pairs of retarded X-operators (Eig. [S|). If the time 
points corresponding to these retarded pairs are (RjT') 
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2T 



2T 



E. Simplified model 


To avoid the sign problem originating from configu¬ 
rations with an odd number of retarded spin-flips/pair- 
hoppings, it is useful to consider a simplified action which 
has a retarded density-density interaction, but only in¬ 
stantaneous spin-flip and pair-hopping terms: 


QSimp 

‘^int 


f dT'na.^(r)(^Ys)f;^'(T-r')n;3,^/(r') 

CH ,(T 0 

13,a' 



dT Jscr 


S{t) + S^t) + P{t) -I- P^(t) 


(37) 


Note that we choose here the screened Hund coupling 
JscT = J{u} = 0) in front of the second term to correctly 
capture the limit of high screening frequency. The re¬ 
tarded density-density interaction of the simplified model 
is 


FIG. 5: Swap updates. We can generate additional configu¬ 
rations with retarded pairs of X operators by swapping the 
end-points of two retarded spin-flip or pair-hopping events. 


and and the configuration before and after the 

swapping is denoted by C and C", respectively, the ac¬ 
ceptance probability for the move is given by 


Pswap(C' C) 


Jr&tiXi Ft) 

JretiT'i — Ft) Jret(Tj ~ Tj) 


(35) 


In particular, this type of update allows us to 
produce configurations with retarded interactions of 
the type -Jret(r' - and - 

We furthermore need the insertion/removal of “kinks” 
-Jretir 'Suppose we have kinks 
of type K and we randomly choose one type for the in¬ 
sertion or removal. For the insertion, we randomly se¬ 
lect the time r for the first operator. If the kink in¬ 
sertion is possible, we compute the length Zmax of the 
interval in which the second operator can be inserted 
and place it at the time t' according to the distribution 
-JretW - t)/ dT'\Jret{T' - t)|. In the reverse 

move, we randomly select one of the + 1 kinks of 

type K, and remove the corresponding X operators if 
there is no other operator in between. The correspond¬ 
ing ratio of acceptance probabilities is 


^ ^kink + 1) — 


«kink + 1 




9,a 

overlap 




dr'l Jret(T' - t)|. 


( 36 ) 




1 — ^bare'^(’^ 

-) + Gret(T), 

(38) 

:-Ar) = 


1 — (^Aiare ~ 

^Jhare')d{T ) 




+ (t^ret('r) 

- 2Jret(r)), 

(39) 


= mlUr] 

1 = (t^bare ~ 

■ 3Tbare)(^('7' ) 




+ (t/ret(T) 

- 3Jret(T)), 

(40) 

3 AT-functions in Eq. 

(05)) become 



KZ3r) 

= Kui-r), 


(41) 

V-Ar) - 

-kZ-At) 

= Ku{t) - 

2KjiT), 

(42) 


= kZA^) 

= Ku(t) - 

iKj{T). 

(43) 


The half-filling condition for the simplified model is 
/^l/2 “ "2 Jscr- 

Because of the spin-dependence of the inter-orbital in¬ 
teraction, Eqs. m and (1321) now read 


Rssi oc e 
Rppt oc c 


-4if j(0-kE'E, 

— 20/Cj(/)+^^ 



•> 


(44) 

(45) 


where I is the distance between the inserted spin-flip or 
pair-hopping operators. Apart from this change, the sim¬ 
ulation proceeds as discussed previously for the case of 
retarded density-density interactions. 


III. TESTS OF THE SOLVER 

We first show the results of some tests of the dou¬ 
ble expansion solver, starting with a model that contains 
only instantaneous spin-flip and pair-hopping terms, and 
a static U. The top panel of Fig. [5] compares the Green’s 
functions obtained with the new solver and with the 
matrix formalism^ for a half-filled orbitally-degenerate 
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FIG. 6: Comparison of simulations results for t/ = 8, J = 8/6, 
^ = 4.5 and /I = 50 (filling n=0.396). The top panel shows the 
Green’s function, with the solid line showing the result from 
the double-expansion solver, while the dashed line has been 
obtained with the matrix formalism. The inset in the top 
panel shows the distribution of hybridization orders, which 
is identical in both methods. The bottom panel shows the 
distribution of spin-flip and, as an inset, the distribution of 
pair hopping orders for this parameter set. 


model with U = 8, J = 1.33, /r = 4.5 and /3 = 50 (DMFT 
solution for a semi-circular density of states with band¬ 
width 4). Both results agree within statistical errors (the 
error bars are comparable to the line thickness), and as 
shown in the inset, also the distribution of the pertur¬ 
bation orders rih is identical. This is because the sam¬ 
pling of the hybridization operators is independent of the 
treatment of the spin-flip and pair-hopping terms (exact 
treatment in the time-evolution in the case of 

the matrix formalism versus stochastic sampling in the 
double-expansion approach). Obviously, also the physi¬ 
cal quantities derived from this distribution of perturba¬ 
tion orders, such as the kinetic energy^^ will agree. 
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FIG. 7: Average perturbation order for Uh, Usf, and riph as a 
function of J for [7 = 3 and half-filling (/3 = 25). 


The bottom panel shows the distribution of spin-flip 
orders Ugi and as an inset the distribution of pair-hopping 
orders Uph- In this simulation, configurations with spin- 
flip orders up to about 25 are relevant, while configu¬ 
rations with more than 3 pair-hopping events are rarely 
generated. The average spin-flip order is 5.93 whereas the 
average pair-hopping order is only 0.18. This is because 
a positive J favors high-spin states (S = 1), and one of 
the high-spin states, + lit)); is an eigenstate of 

(S + S^). For negative J, the situation is opposite and 
the average perturbation order for the pair-hopping term 
becomes large. 

The scaling of the perturbation orders with J is il¬ 
lustrated in Fig. [7] for a half-filled metallic system. For 
small I J|, Usf and riph grows roughly proportional to J^. 
This result is expected since the J-terms are treated by 
a weak-coupling expansion, and in the simulation with 
static J, each S or P operator must be balanced by a 
or pl. Hence, the weight (fTTI) is proportional to (J^)”=f 
and (J^)”p‘*. At larger | J|, deviations from the quadratic 
behavior appear, due to changes in the hybridization or¬ 
der and interference between spin-flip and pair-hopping 
terms. We also notice that Ugf grows more rapidly than 
riph on the J > 0 side, while it is the opposite on the 
J < 0 side. As mentioned above, this is because of the in¬ 
creased weight of high-spin (low-spin) states in the model 
with J > 0 (J < 0). 

We next test the implementation with retarded inter¬ 
actions. As was already mentioned, in the case of a re¬ 
tarded J{t), the expansion in the retarded spin-flip and 
pair-hopping terms produces a sign problem. Therefore, 
we first consider the simplified model ()37l) where only 
the density-density interaction is retarded, whereas the 
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spin-flip and pair-hopping terms are instantaneous. For 
this test, we assume a static interaction U, so that the 
retarded density-density part arises from the retardation 
of J(t). The same-spin density-density interaction thus 
becomes U — 3>7(r), implying = —3Kj{t) for 

a ^ f 3 , and the half-filling condition is ^1/2 = §17—| Jscr- 

For the frequency-dependence, we take a single-boson 
model with a screening frequency ojj and a coupling 
strength Aj, which corresponds to 

J(CC) = Jbare+^^, (46) 

2 17 ^ 

and hence Jscr = Aiare —In Fig. [5] we show re¬ 
sults for 17 = 4, /r = /ii/ 2 , 7bare = 1, fixed Jscr = 0.4 
and different values of uij. In the limit of large wj, 
the results should approach those for the frequency in¬ 
dependent parameters U and Jscr (see inset of the bot¬ 
tom panel). Indeed, as seen in the top panel of Fig. |8l 
the Green’s function approaches the result for a static 
Jscr as LOj oo, and the same is true for the imaginary 
part of the self-energy at the lowest Matsubara frequency 
uji = 7r//3 (inset). Two-particle quantities, such as the 
equal-time spin-correlation function (5^) (middle panel) 
or the double occupancy {rufriii) (bottom panel) also ap¬ 
proach the correct values in the limit of a large screening 
frequency. 

Besides the calculations with instantanteous spin-flip 
and pair-hopping terms, we also show the result from the 
simple density-density approximation, which ignores the 
spin-flip and pair-hopping contributions and treats only 
the (dynamical) J^. As is evident from the plots, these 
two approximations can produce quite different results. 
For example, the density-density calculation exhibits a 
transition to an insulating high-spin state at screening 
frequency ujj « 7, whereas the calculation with instanta¬ 
nteous spin-flip and pair-hopping terms yields a metallic 
solution down to wj = I. 

Finally, let us quantify the effect of the retarded spin- 


TABLE I: Average signs and perturbation orders for a simula¬ 
tion with instantaneous spin-flips and retarded XX operator 
pairs, 17 = 4, Jbare = 1 and Jscr = 0. n-inst is the average num¬ 
ber of spin-flip (S or 5"^) operators, while riret is the average 
number oi XX operators linked by Jret('r). rikink corresponds 
to the subset of kink-operators 


toj 

P 

sign 

^inst 

Tlret 

^kink 

2 

1 

0.6768 

0.379 

1.106 

0.451 

2 

2 

0.2807 

1.139 

2.805 

0.821 

2 

3 

0.1036 

1.913 

4.472 

1.108 

2 

4 

0.0366 

2.665 

6.015 

1.326 

2 

5 

0.0127 

3.397 

7.455 

1.500 

10 

1 

0.6369 

0.432 

1.263 

0.501 

10 

2 

0.2093 

1.388 

3.252 

0.918 

10 

3 

0.0623 

2.363 

5.062 

1.206 

10 

4 

0.0196 

3.290 

6.760 

1.421 

10 

5 

0.0064 

4.175 

8.426 

1.601 




FIG. 8: Green’s functions and local observables for U — 4, 
Jbare = 1, Jscr = 0.4, and different screening frequencies cuj 
{fi = /ii /2 and j3 = 25). Results for the static limit oij —>■ oo 
{U = 4, J = Jscr = 0.4) are also shown (horizontal bars in the 
bottom two panels) or subtracted (top panel). The curves la¬ 
beled by J are from simulations which include instantaneous 
spin-flip and pair-hopping terms (with J = Jscr), while the 
curves labeled by are for the density-density approxima¬ 
tion. Insets: ImE at the lowest Matsubara frequency (top 
panel), and J(tt)) on the real-frequency axis (bottom panel). 
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FIG. 9: Average perturbation order for n^, risf, and riph as 
a function of filling (per spin and orbital) for U — 5,8 and 
J = U/6 iP = 50). 


FIG. 10: Top two panels: Average perturbation orders rih, 
riaf, and riph as a function of the crystal field splitting A for 
U — 5 and J = U/6 {P = 50). The lowest panel shows the 
orbital polarization ni — n 2 . 


IV. RESULTS 


flips and kinks. To simplify the calculations, we switch off 
the pair-hopping terms and show the results after one it¬ 
eration starting from a metallic solution (noninteracting 
hybridization function for P = 50). We consider a model 
with {7 = 4, Jbare = 1 and Jscr = 0. Table U shows the 
average order of instantaneous and retarded XX pairs, 
riinst and riret, as well as the average sign, for screening 
frequencies uij = 2 and loj = 10 and different inverse 
temperatures p. An instantaneous XX pair is either a 
S or 5"^ operator, while Jret can connect different types 
of X operators. The average number of kink-operators 
(a subset of the retarded operators) is listed as 
Ukink- The sign drops exponentially with P and, at least 
in the temperature range considered, faster than expo¬ 
nentially with the average order of retarded XX pairs. 
This limits the simulations of the full model to high tem¬ 
peratures. In the following, we will therefore focus on the 
simplified model, with only instantaneous spin-flips and 
pair-hoppings. 


A. Two-orbital model with static interactions 

We first consider the model with static interactions 
U = 5, J = 5/6 and U = 8, J = 8/6, respectively, and 
plot the average perturbation orders for the hybridiza¬ 
tion, spin-flip and pair-hopping terms as a function of 
filling (per orbital and spin) in Fig. [HI Near half-filling, 
and for the stronger interaction also near quarter-filling, 
the average hybridization order, and hence the kinetic 
energy, is suppressed. This suppression appears to be 
related to the spin-freezing phenomenoni^ As was shown 
in previous studies^i^ also the two-orbital model (both 
with and without spin-flip and pair-hopping terms) ex¬ 
hibits a bad metallic state with “frozen” local moments 
in a certain filling range close to half-filling and quar¬ 
ter filling. The downturns near filling j and ^ in the 
hybridization order coincide with the onset of this spin¬ 
freezing regime (see Fig. S3 in Ref. [2^ and Fig. 12 in 
Ref. El . 
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While the pair-hopping order shows a similar down¬ 
turn, which is also explained by the appearance of S' = 
1 moments, the spin-flip order increases systematically 
with filling and takes its largest value in the half-filled 
Mott insulator. This is because for J > 0, the Mott 
insulator is dominated by high-spin states with one elec¬ 
tron in each orbital. While the spin-flip operators can 
act on some of these spin-triplet states, the pair-hopping 
operator cannot. The situation would be opposite for 
J < 0. 

If the crystal field splitting A is increased in the half- 
filled high-spin Mott insulator, one either observes a first 
order transition to the low-spin insulator (at large U), 
or first a transition into a strongly correlated metallic 
state, followed by a second transition to the low-spin 
insulator.— For the latter case (C/ = 5, J = t//6, /3 = 50) 
we plot the different perturbation orders and the orbital 
polarization as a function of A in Fig. 1101 The hybridiza¬ 
tion order rih increases after the transition into the metal, 
since it is proportional to the kinetic energy, and be¬ 
comes very low in the orbitally polarized low-spin insula¬ 
tor. The spin-flip order risf remains approximately inde¬ 
pendent of A in the high-spin insulator, decreases contin¬ 
uously with increasing A in the metal, and then drops to 
very small values in the low-spin insulating phase. The 
average pair hopping order riph, which is very low in the 
high-spin insulator, grows with increasing orbital polar¬ 
ization in the metal phase, and then jumps to a maximum 
value on the insulating side of the metal-low-spin insula¬ 
tor phase boundary. While riph is large in the low-spin 
insulating phase, it decreases with increasing A, because 
the population of the higher orbital with two electrons 
becomes increasingly costly. 


B. Two-orbital model with dynamically screened J 

In this section, we present results for a two orbital 
model with static U and dynamically screened J. To 
avoid the sign problem associated with retarded spin-flip 
or pair-hopping terms, we only treat the Jz component 
dynamically [Jz{uj) = and approximate the spin- 

flip and pair-hopping terms with a static Jscr = J(y> = 
0). For the frequency dependence we adopt the single bo¬ 
son model (l46l) . Since we will be interested in particular 
also in negative Jscr, our calculations may be viewed as 
simple model calculations for the alkali-doped fullerides. 
In these materials, the effectively negative Jscr favors a 
low-spin state containing an intraorbital electron pair, 
which has been argued to drive the s-wave superconduc¬ 
tivity next to the Mott insulating phase.— Here, we 
will not study superconductivity (for a brief discussion 
of technical aspects related to simulations in the super¬ 
conducting phase, see Appendix but the spin state 
transitions which occur as Jscr is varied. 

The top panel in Fig. [TT] plots (S'^) for a half-filled 
model with [/ = 4, unscreened Jbare = 1, screening fre¬ 
quency wj = 2 and several values of Jscr in the range 



FIG. 11: First two panels: {S\) as a function of Jscr for 1/ = 4, 
Jbare = 1, Wj = 2 (first panel) and loj = oo (second panel), 
M = Mi /2 and J = 25. Both results for the simplified model 
with static spin-flip and pair-hopping terms (“J”) and for the 
density density approximation (“Jz”) are shown. Third and 
fourth panel: average perturbation orders for the spin-flip 
and pair-hopping terms in the simulation with cjj = 2 (third 
panel) and ujj = oo (fourth panel). 


—0.2 < Jscr < 0.5. The curve labeled by ‘J’ shows 
the results of simulations with (static) spin-flip and pair¬ 
hopping terms. For sufficiently large Jscr, the half-filled 
system is in a high-spin insulating state and {SD « 1. 
Around Jscr = 0.43, a transition to a metallic phase 
occurs, and this phase is stable down to Jscr ~ —0.13. 
Only for screened values below —0.13 do we find a low- 
spin insulating phase. The curve labeled by ‘Jj’ shows 
the results obtained from the density-density approxi¬ 
mation (no spin-flip and pair-hopping terms). Consis¬ 
tent with Fig. [51 the stability range of the high-spin in¬ 
sulator is enhanced in the Jz approximation. On the 
other hand, the transition to the low-spin insulator oc¬ 
curs at almost the same Jscr as in the calculation with 
spin-flips/pair-hoppings. This is because in the range 
—0.2 < Jscr 0.2 the average perturbation order for 
pair-hoppings and spin-flips is very low, as shown in the 
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FIG. 12: (Sf) as a function of Jscr for (7 = 4, Jbare = 1, 
ojj = 2, ^ = /ii/ 2 , /3 = 25 and different screening frequencies 


third panel. Note that in the simplified model, Jscr is 
used for the spin-flip and pair-hopping terms, while the 
density-density part is dynamical (i.e. approaches Jbare 
at high frequencies). This explains why for small neg¬ 
ative Jscr, the spin-flip order is slightly larger than the 
pair-hopping order. It also explains the drop in the per¬ 
turbation order across the transition from the metal to 
the high-spin insulator: The simplified model breaks the 
spin SU(2) symmetry and hence the three-fold degener¬ 
acy of the high-spin states. In the high-spin insulating 
phase, the weights of the states |tt) £md |4.|) dominate. 
Consequently, across the transition into the high-spin in¬ 
sulating state, the weights of the states |ti) and Hf) de¬ 
crease. Since the spin-flip operators act on these states, 
the spin-flip perturbation order decreases as one enters 
into the high-spin insulating phase. 

For comparison, we also show the results for the static 
limit (wj = oo) in the second and fourth panel. As ex¬ 
pected, the high-spin insulator becomes less stable, since 
the transition occurs at Jgcr < Jbare- The transition to 
the low-spin insulator exhibits a rather large difference in 
the critical Jscr between the ‘ J’ and ‘ J^’ approximations. 
This can be ascribed to the fact that the average pertur¬ 
bation orders for the pair-hopping and spin-flip interac¬ 
tions becomes larger compared to the wj = 2 case. In 
the Jz approximation, the metallic solution is stabilized 
in the vicinity of the transition to the low-spin insula¬ 
tor. This result is qualitatively similar to the behavior 
of the Holstein-Hubbard model near the transition to the 
bipolaronic insulator, where a large screening frequency 
stabilizes the metallic phase (see e. g. Fig. 2 in Ref. [25h . 
One can explain this correlation effect by translating the 
frequency-dependent density-density interaction into an 
effective bandwidth reduction4^ 

In the calculation with spin-flips and pair-hoppings, 
the effect of wj on the critical Jgcr for the transition to 



FIG. 13: Phase diagram in the space of Jscr and ujj for (7 = 4, 
Jbare = 1, p = /ri /2 uud /3 = 25. The rectangle at tuj = 1 
indicates the parameter values for “fullerene compounds” (see 
text). 

the low-spin insulator is much smaller (the transition oc¬ 
curs at almost the same Jscr for luj = oo and ujj = 2). 
This seemingly small effect is due to a compensation be¬ 
tween the above-described stabilization mechanism for 
the metallic solution and the increase of the pair-hopping 
perturbation orders, which favor the insulating solution. 

For the calculation with spin-flip/pair-hopping terms, 
we plot {S'^) for different screening frequencies ujj in 
Fig. HH With decreasing ujj (< 4), the transition to the 
low-spin insulator shifts to less negative values of Jscr, 
while the critical value for the transition to the high-spin 
insulator shows a non-uniform behavior. Also, in the 
limit of small wj, the transition to the low-spin insulator 
is marked by a large jump in {S'^). The dependence of the 
phase boundaries on uij is plotted over a wider range of 
uij in Fig. [131 As ujj is reduced from the large-frequency 
limit, both phase boundaries shift to smaller Jscr, be¬ 
cause the increasing effect of the Jbare = 1 on fast spin 
fluctuations stabilizes (destabilizes) the high-spin (low- 
spin) insulator. At smaller wj, the band renormalization 
effect, which enhances the stability of the low-spin insu¬ 
lator, reverses this trend. 

While a discussion of fullerene compounds based on 
results for a two-orbital model is dangerous, because the 
fullerides are three-orbital systems, it is still interesting 
to comment on the realistic values of wj, Jscr, and (7, 
if translated to the current set-up (with bandwidth 4): 
(jjj ^ 1, Jscr between —0.15 and —0.1, and U ^ 5-10i^2i2S 
Our results in Fig. (T^j and [T3| suggest that the alkali- 
doped fullerides are located near the transition between 
the metal and low-spin insulator (see box in Fig. I13p . in 
agreement with experiments4^“— 

The effect of a larger U on the phase diagram of Fig.[T51 
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is to enhance the stability region of the high-spin insula¬ 
tor, and to a lesser extent also of the low-spin insulator. 
For example, for U = 6 , the transitions to the high-spin 
(low-spin) insulator occur at Jscr = 0.13 (Jscr = —0.07) 
for LUj = oo, and at Jscr = 0.06 (Jscr = —0.03) for wj = 2. 


V. CONCLUSIONS AND OUTLOOK 

We have presented a double-expansion algorithm for 
multi-orbital impurity models which combines a hy¬ 
bridization expansion with a weak-coupling expansion 
for the spin-flip and pair-hopping terms. This algo¬ 
rithm is based on the economical segment representation 
and avoids computationally expensive matrix multiplica¬ 
tions. By construction, it performs particularly well in 
the regime of weak Hund coupling. 

A main motivation for introducing this algorithm is 
the ability to treat dynamically screened J(w). We have 
explained the Monte Carlo procedures for the sampling 
of retarded density-density, spin-flip and pair-hopping 
terms. In practice, however, the sign problem associ¬ 
ated with retarded spin-flips and pair-hoppings forced 
us to consider a simplified model, in which only the Jz 
component is treated as dynamical, while spin-flip and 
pair-hopping terms are approximated as static operators. 
This approximation allows sign-free simulations of the 
two-orbital model, and to explore the effect of a dynam¬ 
ically screened J{uj) in a set-up which recovers the rota- 
tionally invariant interaction in the limit of large screen¬ 
ing frequency. 

Extending our formalism to models with more than 
two orbitals requires additional Monte Carlo moves which 
change the “winding number” of configurations (Ap¬ 
pendix]^. These updates will introduce negative weight 
configurations even in the case of the simplified model. 
Keeping track of the complex topology of the segment 
configurations for models with more than two orbitals 
becomes challenging. It may thus be simpler to pur¬ 
sue an alternative approach, which is not as efficient 
as the segment-based algorithm, but more flexible and 
easier to implement: As explained in Appendix |B] one 
can rewrite the spin-flip and pair-hopping operators us¬ 
ing auxiliary spin variables. In combination with the 
hybridization expansion algorithm in the matrix formu¬ 
lation and a weak-coupling expansion in the spin-flip and 
pair-hopping terms, this auxiliary field sampling will al¬ 
low to generate configurations with an odd number of 
operators, and hence non-zero winding number. How¬ 
ever, the physical configurations are in this approach al¬ 
ways combined with fictitious ones, whose contributions 
to the Monte Carlo measurements average to zero, but 
exacerbate the sign problem. 

On the application side we have presented results for 
the spin state transitions in a two-orbital model with 
static U and dynamically screened J(w). These results 
show that our simplified model with static spin-flip and 
pair-hopping terms produces results which differ signifi¬ 


cantly from the density-density approximation near the 
transition to the high-spin Mott insulating state. On the 
other hand, the transition to the low-spin insulating state 
occurs at small negative Jscr- Because of this, the pertur¬ 
bation orders for the spin-flip and pair-hopping terms are 
small, which results in almost identical transition points 
as in the density-density case. At least in our single¬ 
boson model, the low-spin insulating phase is stabilized 
significantly in the limit of low screening frequency ujj. 
Choosing realistic parameter values for fullerides, we find 
that within our two-orbital description, the system is on 
the verge of the transition from the metal to the low-spin 
insulating phase. 

Because of the relevance of multi-orbital models with 
overscreened J for the physics of alkali-doped fullerides, 
an interesting future application will be the study of 
superconductivity in this model. In particular, it will 
be possible to investigate the effect of the frequency- 
dependent Hund’s coupling on the properties of the su¬ 
perconducting state. To map out the stability regions of 
superconducting or other symmetry-broken phases, it is 
sufficient to measure appropriate susceptibilities in the 
symmetric state, using the algorithm discussed in this 
paper. In order to enter the superconducting phase, ad¬ 
ditional Monte Carlo updates are needed. We briefly 
discuss this issue in Appendix [C] Testing the efficiency 
of the proposed method in the superconducting state is 
an important and interesting future problem. 
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Appendix A: Three-orbital model 

For multi-orbital models with more than two orbitals, 
the updates described in Sec. Ill Cl are not sufficient for 
an ergodic sampling. This is because configurations with 
non-zero “winding number” contribute to the partition 
function, and these configurations cannot be generated 
by the insertion/removal of operator pairs SS^, S or 
PP\ P^P. We will focus the discussion here on the 
three-orbital case and the spin flip operators; the gener¬ 
alization to more orbitals and to the pair-hopping case is 
straightforward. 

In the three orbital context, it is natural to work with 
the spin-flip operators Sap = /3 = 1) 2, 3 

and a ^ (3). The number of these operators will be 
denoted by ns- An example of a ns = 3 configuration 
with operators 8 ^ 2 , S 21 and S '13 is illustrated in the top 
panel of Fig. [TH 




Appendix B: An alternative way to perform the 
double expansion: “auxiliary spin” method 
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Here, we propose an alternative method to treat spin- 
flip and pair-hopping interactions. For simplicity, let us 
consider static spin-flip and pair-hopping terms. In the 
scheme discussed in Sec. Ill Bl for the two-orbital model 
with static spin-flip and pair-hopping interactions, the 
corresponding perturbation orders are always even: The 
S operator is always paired with a operator and the 
same is true for P and . 

Alternatively, we can perform the simulation by in¬ 
troducing auxiliary spins, similarly to the strategy em¬ 
ployed in Ref. for the interaction-expansion impurity 
solveri^ In this alternative method, we rewrite an O op¬ 
erator {O = S, ,P, and P^) as 

E (Bl) 

Sa=±l 


11 


0 


where 


FIG. 14; Top panel: Configuration with ns = 3 spin flips and 
nonzero winding number in the three-orbital model. Bottom 
panel: Configuration with ns = 2 spin-flips from which the 
configuration in the upper panel can be obtained using the 
procedure described in Appendix m 


To produce this configuration, we can start from a 
neighboring 532, 5'23 pair of spin flips as shown in the bot¬ 
tom panel. These operators define an interval of length 
^max in which we randomly choose the time r. We then 
propose to insert a spin flip S'21 at r, and simultaneously 
replace the S23 operator by an S'13 operator. The pro¬ 
posal probability for this move is p^’^°^{ns ns + 1) = 
(Of course, the configuration in orbital 1 must 
be compatible with this operator insertion and operator 
replacement, otherwise the move is rejected.) For the in¬ 
verse move, we randomly select a spin-flip operator and 
propose to remove it by simultaneously changing the type 
of the spin-flip operator to the right. In this case, the pro¬ 
posal probability is p^’^°^{ns -I- 1 —S' ns), and the ratio of 
acceptance probabilities becomes 


i?(ns^ns + l) = -^^^%^. (Al) 

ns + 1 


Similarly, we could have proposed to insert a S'13 opera¬ 
tor and to replace the S23 by an S21 (which is possible 
only if the If state is occupied in orbital 1). Note the 
minus sign in the acceptance ratio, which indicates that 
the sampling of the configurations with nonzero winding 
numbers introduces a sign problem. 


ds ^=0 + Sail (B2) 

with Sa the auxiliary spin, 7 a positive real number, and 
/ the identity operator. We then expand the partition 
function in powers of the hybridization operators and the 
operators and perform the sum over the spin ori¬ 
entations Sa = if in the Monte Carlo sampling. 

In this scheme, in addition to the insertion/removal 
of pairs of hybridization operators, we use the following 
updates: 

1. insertion and removal of ^Ss^{ts), 

2. insertion and removal of 

3. insertion and removal of ^Fs^(Tp), 

4. insertion and removal of ^Pj^{Tp)- 

Note that, in this method, the perturbation order can be 
odd because of the Sall term. In evaluating the config¬ 
uration weights, if we employ the segment formalism, we 
have to sum up the weights for 2^ configurations with 
N being the total perturbation order of the ^Osa op¬ 
erators. This is because at each vertex, we have 
the choice between O and Sail operators. It is easy 
to show that many of these 2^ configurations have zero 
weight. Therefore, in practice, we can reduce the number 
of configurations in the calculation of the weight. An¬ 
other possibility is to employ the matrix formalismi^ In 
this case, each of the operators can be represented 
by a single matrix, and the summation of the weights 
over 2 ^ configurations can be avoided. In practice, if we 
use a block-diagonalized matrix representation and lo¬ 
cal conserved quantities we only need to apply the 
Sail operator to block-diagonalized matrices which con¬ 
tain nonzero matrix elements of the operator O. A larger 
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FIG. 15: Illustration of a configuration with one pair-hopping 
operator and two anomalous hybridization functions A^'*' 
which contributes to the partition function of a two-orbital 
model with intra-orbital singlet pairing. This configuration 
can be obtained from the one with two full lines in orbital 
a? = 1 by the procedure described in Appendix El 

7 improves the acceptance ratio of the above-mentioned 
updates; however, a larger 7 also produces a larger num¬ 
ber of samples with negative weight. Therefore, one has 
to find the optimal value for 7 , which will depend on e.g. 
the size of Jscr- 

In the presence of retarded density-density interac¬ 
tions, the configuration weight also contains a bosonic 
factor (Sec. Ill D II) . In the auxiliary spin algorithm, we 
evaluate this bosonic factor with the operators O instead 
of . This procedure produces unphysical weights for 
configurations involving Sa'yl operators due to an incon¬ 
sistency in evaluating the local trace and the bosonic fac¬ 
tor. However, these contributions will average to zero in 


the Monte Carlo sampling, and thus we still get correct 
thermal averages for the physical quantities. Because of 
the “fictitious” weights involving SajI operators, it is 
clear that the auxiliary spin method will have a more 
severe sign problem, and thus be less efficient than the 
scheme described in Sec. Ill Cl However, it circumvents 
the algorithmic complexity originating from the increase 
in the number of orbitals, as described in Appendix El 


Appendix C: Simulations in the superconducting 
state 

In a simulation of a symmetry-broken phase, we some¬ 
times need special updates in addition to the updates 
used in the simulation of the normal phase, as was 
pointed out in Ref. For example, let us consider 
the case of s-wave superconductivity in a two-orbital 
model with a negative static Jscr- A negative Jscr fa¬ 
vors intraorbital singlet pairing4i“— To sample this 
superconducting state, one will need additional updates 
such as the insertion and removal of a P [P^] operator 
and two anomalous hybridization functions (see Fig. 1151) . 
These anomalous hybridization functions correspond to 
two Hhyh operators with flavor (l,t) and (l,-() [( 2 ,t) 
and (2,4.)] and two operators with flavor (2,'|') and 
(2,4.) [(l,t) and (l,i)]. Alternatively, one could think 
of an update in which one instantaneous P [P^] opera¬ 
tor is replaced by two Hhyh operators with flavor ( 2 ,'|') 
and (2,4,) [(l,t) and (1,4-)] and two P^yt operators with 
flavor (l,t) and (l,s() [(2,t) and (2,4,)]. 
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